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. We introduce and investigate the stochastic dynamics of the density of local extrema (minima 

and maxima) of non-equilibrium surface fluctuations. We give a number of exact, analytic results 
for interface fluctuations described by linear Langevin equations, and for on-lattice, solid-on-solid 

■ surface growth models. We show that in spite of the non-universal character of the quantities studied, 
their behavior against the variation of the microscopic length scales can present generic features, 
characteristic to the macroscopic observables of the system. The quantities investigated here present 
us with tools that give an entirely un-orthodox approach to the dynamics of surface morphologies: 
a statistical analysis from the short wavelength end of the Fourier decomposition spectrum. In 

■ addition to surface growth applications, our results can be used to solve the asymptotic scalability 
O ' problem of massively parallel algorithms for discrete event simulations, which are extensively used 

in Monte- Carlo type simulations on parallel architectures. 



C3 . I. INTRODUCTION AND MOTIVATION 
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The aim of statistical mechanics is to relate the macroscopic observables to the microscopic properties of the 
system. Before attempting any such derivation one always has to specify the spectrum of length-scales the analysis 
Ch ■ will comprise: while 'macroscopic' is usually defined in a unique way by the every-day-life length scale, the 'microscopic' 
' is never so obvious, and the choice of the best lower-end scale is highly non-universal, it is system dependent, usually 
(— I , left to our physical 'intuition', or it is set by the limitations of the experimental instrumentation. It is obvious that 
Q ' in order to derive the laws of the gaseous matter we do not need to employ the physics of elementary particles, 
O , it is enough to start from an effective microscopic model (or Hamiltonian) on the level of molecular interactions. 
Then starting from the equations of motion on the microscopic level and using a statistic and probabilistic approach, 
the macroscale physics is derived. In this 'long wavelength' approach most of the microscopic, or short wavelength 
^ information is usually redundant and it is scaled away. 

CO Sometimes however, microscopic quantities are important and directly contribute to macroscopic observables, e.g., 

\l the nearest- neighbor correlations in driven systems determine the current, in model B the mobility, in kinetic Ising 
' ' ' model the domain-wall velocity, in parallel computation the utilization (efficiency) of conservative parallel algorithms, 
] etc. Once a lower length scale is set on which we can define an effective microscopic dynamics, it becomes meaningful 
' to ask questions about local properties at this length scale, e.g. nearest neighbor correlations, contour distributions, 
] extremal-point densities, etc. These quantities are obviously not universal, however their behavior against the variation 
. of the length scales can present qualitative and universal features. Here we study the dynamics of macroscopically 
' rough surfaces via investigating an intriguing miscroscopic quantity: the density of extrema (minima) and its finite 
size effects. We derive a number of analytical results about these quantities for a large class of non-equilibrium surface 
fluctuations described by linear Langevin equations, and solid-on-solid (SOS) lattice-growth models. Besides their 
' ^ obvious relevance to surface physics our technique can be used to show [Q the asymptotic scalability of conservative 
massively parallel algorithms for discrete-event simulation, i.e., the fact that the efficiency of such computational 
schemes does not vanish with increasing the number of processing elements, but it has a lower non-trivial bound. 
The solution of this problem is not only of practical importance from the point of view of parallel computing, but it 
has important consequences for our understanding of systems with asynchronous parallel dynamics, in general. There 
are numerous dynamical systems both man-made, and found in the nature, that contain a "substantial amount" of 
}^ ' parallelism. For example, 
5^ . 1) in wireless cellular communications the call arrivals and departures are happening in continuos time (Poisson 
arrivals), and the discrete events (call arrivals) are not synchronized by a global clock. Nevertheless, calls initiated in 
cells substantially far from each other can be processed simultaneously by the parallel simulator without changing the 
poissonian nature of the underlying process. The problem of designing efficient dynamic channel allocation schemes 
for wireless networks is a very difficult one and currently it is done by modelling the network as a system of interacting 
continuous time stochastic automata on parallel architectures 

2) in magnetic systems the discrete events are the spin flip attempts (e.g., Glauber dynamics for Ising systems). 
While traditional single spin flip dynamics may seem inherently serial, systems with short range interactions can 
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be simulated in parallel: spins far from each other with different local simulated times can be updated simultane- 
ously. Fast and efficient parallel Monte-Carlo algorithms are extremely welcome when studying metastable decay and 
hysteresis of kinetic Ising ferromagnets below their critical temperature, see |^ and references therein. 

3) financial markets, and especially the stock market is an extremely dynamic, high connectivity network of relations, 
thousands of trades are being made asynchronously every minute. 

4) the brain. The human brain, in spite of its low weight of approx. 1kg, and volume of 1400 cm"^, it contains 
about 100 billion neurons, each neuron being connected through synapses to approximately 10,000 other neurons. 
The total number of synapses in a human brain is about 1000 trillion (10^^). The neurons of a single human brain 
placed end-to-end would make a "string" of an enourmous lenght: 250,000 miles Assuming that each neuron of a 
single human cortex can be in two states only (resting or acting), the total number of different brain configurations 
would be 2^" . According to Carl Sagan, this number is greater than the total number of protons and electrons of 
the known universe, Q. The brain does an incredible amount of parallel computation: it simultaneously manages all 
of our body functions, we can talk and walk at the same time, etc. 

5) evolution of networks such as the internet, has parallel dynamics: the local network connectivity changes concur- 
rently as many sites are attached or removed in different locations. As a matter of fact the physics of such dynamic 
networks is a currently heavily investigated and a rapidly emerging field |^ . 

In order to present the basic ideas and notions in the simplest way, in the following we will restrict ourselves to one 
dimensional interfaces that have no overhangs. The restriction on the overhangs may actually be lifted with a proper 
paramctrization of the surface, a problem to which we will return briefly in the concluding section. The first visual 
impression when we look at a rough surface h{x, t) is the extent of the fluctuations perpendicular to the substrate, in 
other words, the width of the fluctuations. The width (or the rms of the height h of the fluctuations) is probably the 
most extensively studied quantity in interface physics, due to the fact that its definition is simply quantifiable and 
therefore measurable: 



w{L,t) - [h{x,t)Y- [h{x,t^ , (1) 

where the overline denotes the average over the substrate. It is well-known that this quantity characterizes the 
long wavelength behavior of the fluctuations, the high frequency components being averag ed out in (§). The short 
wavelength end of the spectrum has been ignored in the literature mainly because of its non-universal character, and 
also because it seemed to lack such a simple quantifiable definition as the width w. 

In the following we will present a quantity that is almost as simple and intuitive as the width w but it characterizes 
the high frequency components of the fluctuations and it is simply quantifiable. For illustrative purposes let us 
consider the classic Weierstrass function defined as the M oo limit of the smooth functions W^Af(a, b; x): 



M 

Wia,b;x)= lim WM{a,b;x)^ lim V a"" cos (6"a;) , a,b>l (2) 



m=0 



Figure la shows the graph of W]\j at a = 2 and 6 = 3 (arbitrary values) for M = 0, 1, 2, 3, 4 in the interval x G [0, 47r]. 
As one can see, by increasing AI we are adding more and more detail to the graph of the function on finer and finer 
length-scales. Thus M plays the role of a regulator for the microscopic cut-off length which is 6^*^, and for M = oo 
and b > a, the function becomes nowhere differentiable as it was shown by Hardy Q . 

Comparing the graphs of Wm for lower M values with those for higher Ms we observe that the width effectively 
does not change, however the curves look qualitatively very different. This is obvious from (^: adding an extra 
term will not change the long wavelength modes, but adds a higher frequency component to the Fourier spectrum 
of the graph. We need to operationally define a quantifiable parameter which makes a distinction between a much 
'fuzzier' graph, such as Af = 4 and a smoother one, such as Af = 1. The natural choice based on Figure la is the 
number of local minima (or extrema) in the graph of function. In Figure lb we present the number of local minima 
um vs. M for two different values oi b, b = 2.8 and b — 1.8, while keeping a at the same value of a = 2. For all b 
values (not only for these two) the leading behaviour is exponential: um ^ A*^. The inset in Figure lb shows the 
dependence of the rate A as a function of b for fixed a. We observe that for 6 > a, A = 6, but below b ~ a the 
dependence crosses over to another, seemingly linear function. For b > a the amplitude of the extra added term in 
Wm+1 is large enough to prevent the cancellation of the newly appearing minima by the drop in the local slope of 
Wm- At b < a the number of cancelled extrema starts to increase drastically with an exponential trend, leading to 
the crossover seen in the inset of Figure 2b. It has been shown that the fractal dimension of the Weierstrass function 
for 6 > a is given by Dq — In 5/ In a |^], ||]. For b < a the Weierstrass curve becomes non- fractal with a dimension 
of Do — 1- By varying b with respect to a, we are crossing a fractal-smooth transition at 6 = a. The very intriguing 
observation we just come across is that even though we are in the smooth regime (6 < a) the density of minima is 
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still a diverging quantity (the b = 1.8 curve in Figure lb). It is thus possible to have an infinite number of 'wrinkles' 
in the Weierstrass function without having a diverging length, without having a fractal in the classical sense. The 
transition from fractal to smooth, as b is lowered appears as a non-analyticity in the divergence rate of the curve's 
wrinkledness. A rigorous analytic treatment of this problem seems to be highly non-trivial and we propose it as an 
open question to the interested reader. 
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FIG. 1. a) the function WM{a, b; x) at a = 2, b = 3 and M = 0, 1, 2, 3, 4. b) The scaling um ~ A*^ of the number of local 
minima for b — 2.8 (plusses) and b = 1.8 (crosses). The inset shows A vs. b while keeping a constant, a = 2. 

The simple example shown above suggests that there is novel and non-trivial physics lying behind the analysis of 
extremal point densities, and it gives extra information on the morphology of interfaces. Given an interface h{x, t), we 
propose a quantitative form that characterizes the density of minima via a 'partition-function' like expression, which 
is hardly more complex than Eq. (^ and gives an alternate description of the surface morphology: 

Uq{L, [h]) = — [K{xi)]'' , q > , Xi are non-degenerate minima of h (3) 

i 

with K{xi) denoting the curvature of h at the local (non-degenerate) minimum point Xi. The variable q can be 
conceived as 'inverse temperature'. Obviously, for q — Q we obtain the number of local minima per unit interface 
length. The rigorous mathematical description and definitions lying at the basis of (|^) is being presented in Section 
IV. The quantity in Eq. (^) is reminiscent to the partition function used in the definition of the thermodynamical 
formalism of one dimensional chaotic maps Q and also to the definition of the dynamical or Renyi entropies of these 
chaotic maps. In that case, however the curvatures at the minima are replaced by cylinder intervals and/or the visiting 
probabilities of these cylinders. 

We present a detailed analysis for the above quantity in case of a large class of linear Langevin equations of type 
dh/dt = -i/(-V2ft)^/2 T^{x,t), where 77 is a Gaussian noise term, and z a positive real number. These Langevin 
equations are found to describe faithfully the fiuctuations of monoatomic steps on various substrates, see for a review 
Ref. jl^ . One of the interesting conclusions we came to by studying the extremal-point densities on such equations is 
that depending on the value of z the typical surface morphology can be fractal, or locally smooth, and the two regimes 
are separated by a critical z value, Zc- In the fractal case, the interface will have infinitely many minima and cusps 
just as in the case of the nondiffercntiable Weierstrass function (||), and the extremal point densities become infinite, 
or if the problem is discretized onto a lattice with spacing a, a power-law diverging behavior is observed as a ^ 0"*" 
for these densities. This sudden change of the 'intrinsic roughness' may be conceived as a phase transition even in an 
experimental situation. Changing a parameter, such as the temperature, the law describing the fiuctuations can change 
since the mechanism responsible for the fluctuations can change character as the temperature varies. For example, 
it has been recently shown using Scanning Tunneling Microscope (STM) measurements [Q, that the fiuctuations of 
single atom layer steps on Cu (111) below T = 300 °C correspond to the perifery diffusion mechanism (z = 4), but 
above this temperature (such as T = 500 °C in their measurements) the mechanism is attachement-detachment where 
z = 2, see also Ref. |l|,|l|. 

The paper is organized as follows: In Sections II and III we define and investigate on several well known on-lattice 
models the minimum point density and derive exact results in the steady-state {t 00) including finite size effects. As 
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a practical application of these on-lattice results, we briefly present in Subsection III.B a lattice surface-growth model 
which exactly describes the evolution of the simulated time-horizon for conservative massively parallel schemes in 
parallel computation, and solve a long-standing asymptotic scalability question for these update schemes. In Section 
IV we lay down a more rigorous mathematical treatment for extremal point densities, and stochastic extremal point 
densities on the continuum, with a detailed derivation for a large class of linear Langevin equations (which are in fact 
the continuum counterparts of the discrete ones from Section II). The more rigorous treatment allows for an exact 
analytical evaluation not only in the steady-state, but for all times! We identify novel characteristic dimensions that 
separate regimes with divergent extremal point densities from convergent ones and which give a novel understanding 
of the short wavelength physics behind these kinetic roughening processes. 



II. LINEAR SURFACE GROWTH MODELS ON THE LATTICE 



In the present Section we focus on discrete, one dimensional models from the linear theory of kinetic roughening 
Let us consider a one dimensional substrate consisting of L lattice sites, with periodic boundary conditions. 
For simplicity the lattice constant is taken as unity, which clearly, represents the lower cut-off length for the effective 
equation of motion. For the moment let us study the discretized counterpart of the general Langevin equation that 
describes the linear theory of Molecular Beam Epitaxy (MBE) ||l^,|6| : 



dMt) = vy^h,{t) - K\/*h,{t) + ri,{t) , (4) 
where rji {t) is Gaussian white noise with 

{'n,it)7^,{t')) =2Dd,^j S{t-t') , (5) 

and is the discrete Laplacian operator, i.e., V^/j = fj+i + fj-i ~ '^fjj applied to an arbitrary lattice function fj. 
This equation arrises in MBE with both surface diffusion mechanism (the 4th order, or curvature term) and desorption 
mechanism (the 2nd order, or diffusive term) present and it has been studied extensively by several authors |]l5| , |l7| . 
Stability requires v > and k > (as a matter of fact, on the lattice is enough to have v > and k > —v/2). Starting 
from a completely flat initial condition, the interface roughens until the correlation length ^ reaches the size of the 
system S, — L, when the roughening saturates over into a steady-state regime. The process of kinetic roughening is 
controlled by the intrinsic length scale [0 , -y/ k/v. Below this lengthscale the roughening is dominated by the surface 
diffusion or MuUins ||l^ term (the 4th order operator) but above it is characterized by the evaporation piece (the 
diffusion) or Edwards- Wilkinson |Q term. Since eq. (^ is translationaly invariant and linear in h it can be solved 
via the discrete Fourier-transform: 

^fc = Ee"'''^*' ^ = ^' " = 0,1,2,... (6) 

i=o 

Then eq. (||) translates into 

dthk{t) = - [2v{l ~ cos(fc)) -f 4k(1 - co%(k)f\ hit) + m{t) , (7) 

with 

{Ut)m' it')} - 2DL 5^^^^,^ 2.,o '^(^ - (8) 
Following the definition of the equal-time structure factor for S{k,t), namely 

SHk,t)LS^^^,,^ mod2.,o = Chkit)h'{t)) , (9) 
one obtains for an initially flat interface: 

5''(fc,<) = S'^ik) (l - e-(4^(i-™«W)+8«(i-'=°<fe))')*^ . (10) 

In the above equation 

S'^ik) = lim S''ik,t) = — ,,,, ,^ -— (11) 

^ ^ i^oo ^ ' 2i/(l-cos(fc)) + 4K(l-cos(/c))2 ^ ' 
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is the steady-state structure factor. 

For V ^ and in the asymptotic scahng Umit where L 3> \/ n/v, the model belongs to the EW universality 
class and the roughening exponent is a = 1/2 (it is defined through the scaling L^" of the interface width {w\{t)) — 

0-lL){TA=i[hi{t)-hf) in the steady state). The presence of the curvature term does not change the universal scaling 
properties for the surface width, and one finds the same exponents as for the pure EW (k = 0) case in eq. (^. For 
V = Q the surface is purely curvature driven {z = 4) and the model belongs to a different universality class where the 
steady-state width scales with a roughness exponent of a = 3/2. 

In the following we will be mostly interested in some local steady-state properties of the surface hi. In particular, we 
want to find the density of local minima for the surface described by . The operator which measures this quantity 
is 



1 ^ 

-^e(/i,_i-/iOe(/i.+i-/iO ■ (12) 



1=1 



This expression motivates the introduction of the local slopes, <j)i = /i^+i — hi. In this representation the operator for 
the density of local minima (for the original surface) is 



1 ^ 

-5]e(-0,_i)e(</..) , (13) 



1=1 



and its steady-state average is {u) — (0(— (/)i_i)8(0i)) ~ (8(— (/)i)8(02)), due to translational invariance. The average 
density of local minima is the same as the probability that at a randomly chosen site of the lattice the surface exhibits 
a local minimum. It is governed by the nearest-neighbor two-slope distribution, which is also Gaussian and fully 
determined by {<f)\) ~ {cf)^) and (0i02): 



where 



P""((/)i,(/)2) cx e-i^^-^"?"*'- , j,fc = 1,2 , (14) 



A- = [ ] (15) 



^1/ 



2\ 



As we derive in Appendix A, the density of local minima only depends on the ratio ((/>i(/)2)/('/'i 

(.).i,arccosf%^) . (16) 



2tt 



Finite-size effects of (u) are obviously carried through those of the correlations. First we find the steady-state structure 
factor for the slopes. Since = (1 ^ e^''^)ft,fc, we have S"^(fc) = 2(1 — cos(fc))S''*(fc). Then from (O) one obtains: 



S"^(fc) = — , for /j 7^ , and S"^(fc) = , for fc = . (17) 

V + 2k(1 — cos(k)) 



The former automatically follows from the X^i^i 0i = relation. Then we obtain the slope-slope correlations 



L-l 



C^(0 ^ = \Y. ^'"^'S^ ( ^ ) ■ (18) 



n=l ^ 

With the help of Poisson summation formulas, in Appendix B we show a derivation for the exact spatial correlation 
function, which yields 

C^(/)^^(^-^i + ^4i^|, |/|<X, (19) 

with 



2k , , 1 - Vl - 

and b EE . (20) 
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We have \a\ < 1 and 6 < 1. The second term in the bracket in ( [i9| ) gives a uniform power law correction, while 
the third one gives an exponential correction to the correlation function in the thermodynamic limit. For ^ and 
L ^ oo one obtains 

Ct{l) = 7^ = ^ , (21) 

where we define the correlation length of the slopes for an infinite system as: 

^ -R5) ■ ^''^ 

In the V ^ Q limit it becomes the intrinsic correlation length which diverges as cJ^ \J k/v and 

In this limit the slopes (separated by any finite distance) become highly correlated, and one may start to anticipate 
that the density of local minima will vanish for the original surface {hi\. In the following two subsections we investigate 
the density of local minima and its finite-size effects for the Edwards- Wilkinson and the MuUins cases. 



A. Density of local minima for Edwards- Wilkinson term dominated regime 

To study the finite size effects for the local minimum density, we neglect the exponentially small correction in ([is|), 
so in the asymptotic limit, where L 3> Cf^{l) decays exponentially with uniform finite-size corrections: 



iy + 2K { Vl - a2 1 - a L 
This holds for the special case k = as well, (in fact, there the exponential correction exatly vanishes) leaving 

- L 

Now, emplying eq. (p^), we can obtain the density of minima as 



Ctil) = -{ko-T] ■ (25) 



1 I Cl{l)\ _ 1 , 1 jl-b /1 + al 



{u)l ^ — arccos — r ~ — arccos(&) + —\ \ / — , (26) 

^' 27r \Cl{Q)j 27r ^ ' 2t:\ I + I - a L ' ^ ' 

Again, for the n — case one has a compact exact expression and the corresponding large L behavior: 

1 / 1 A 1 1 1 

<">-=2^^'-^^°%-L^j"4 + 2^Z' (2^) 

which can also be obtained by taking the k ^ limit in (p6|). To summarize, as long a,s v ^ 0, the model belongs to 
the EW universality class, and in the steady state, the density of local minima behaves as 

, ^ , ^ const. 

{u)l ^ {u)oo + , (28) 

where {u)oc is the value of the density of local minima in the thermodynamic limit: 

(m)oo = 77- arccos(fo) . (29) 
Zn 

Note that this quantity can be small, but does not vanish if v is close but not equal to 0. Further, the system exhibits 
the scaling (|28| ) for asymptotically large systems, where L ^ It is important to see in detail how (zi)oo behaves 
as 0: 

^ ™ (l - V2VT^) ^ i- arccos [l _ i- [2 ^ • (^0) 



Thus, the density of local minima for an infinite system vanishes as we approach the purely curvature driven {v — > 0) 
limit. Simply speaking, the local slopes become "infinitely" correlated, such that C'^{1) diverges [according to eq. 
(p3|)], and the ratio C'^(Z)/C^(0) for any fixed / tends to 1. This is the physical picture behind the vanishing density 
of local minima. 
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B. Density of local minima for the Mullins term dominated regime 



Here we take the ly limit first and then study the finite size effects in the purely curvature driven model. The 
slope correlations are finite for finite L as can be seen from eq. (p^, since the n — term is not included in the sum! 
Thus, in the exact closed formula a careful limiting procedure has to be taken which indeed yields the internal 
cancellation of the apparently divergent terms. Then one obtains the exact slope correlations for the 1^ = case: 



and for the local minimum density: 



1 / 6 \ V3 1 

(.), = -arccos(^l--^jc.-— (32) 



It vanishes in the thermodynamic limit, and hence, one observes that the limits v and L ^ oo arc intcrchangable. 

For = 0, i is directly associated with the correlation length and we can define = L/Q. Then the correlations and 
the density of local minima takes the same scaling form as eqs. (|2^) and (^0[): 

ct{i)^§:{e^-\i\) , (33) 

and 

{u),o.^^. (34) 



C. Scaling considerations for higher order equations 

Let us now consider another equation but with a generalized relaxational term that includes the Edwards Wilkinson 
and the noisy Mullins equation as particluar cases: 

dMt) = -iy{-\/Y^\ + 7j,{t) . (35) 

where z is a positive real number (not necessarily integer). Other z values of experimental interest are z — 1, relaxation 
through plastic flow, p9| , p^ ), and z = 3 terrace-diffusion mechanism [pi] . For early times, such that t ^ L^, the 
interface width (wj^ (i)yincreases with time as 

K(^))^^^^ (36) 

where P — (z — l)/2z |p^ , |l8| . In the t ^ oo limit, where t ^ L^, the interface width saturates for a finite system, 
but diverges with L according to {wj^{oo)) ^ L^°' where a = (z — l)/2 is the roughness exponent ^4 18|. 

For z = 4 (curvature driven interface) we saw that the slope fluctuation behaves as C^(0) = ^ L. For higher 
z for the slope-slope correlation function one can deduce 

L-l 2^, 

CtiD-jj: — ^. (37) 

^„tl.[2(l-cos(2fn))]- 

It is divergent in the L — > oo limit, as a result of infinitely small wave- vectors ~ and we can see that 

It is also useful to define the slope difference correlation function 

Gt(/)^((0.+i-0.)') (39) 
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for which one can write 



^ ^lv[2{l~cos{^^))]— 



For the small wave- vector behavior we can again deduce that for z > 5 

Glil) ^ L'-"^ f . (41) 

One may refer to this form as "anomalous" scaling for the slope difference correlation function in the following 
sense. For z < 5 the scaling form for G^(0 follows that of C^(0) [eq. (||)], i.e., 6^(0 - Z^-^. For z > 5 [eq. it 
obviously features a different I dependence and an additional power of L, and it diverges in the L — > oo limit. 

Having these scaling functions for large L, we can easily obtain the scaling behavior for the average density of local 
minima. Exploiting the identity 



C^(/) = C^(0)--Gt(0 (42) 



we use the general form for the local minimum density: 



1 (C'^l{^)\ 1 1 (, const. ^^ 1 



— ■ arccos — ^ = — arccos 1 ^ ~ — arccos 1 ^ ~ — (43) 

27r \Ct{Q)l 2^ \^ 2C7^(0)y 2ti \ \? ) L ^ ' 

Note that this is the scaling behavior for all z > 5. It simply shows the trivial lower bound for {u): since there is 
always at least one minima (and one maxima) among the L sites, it can never be smaller than 1/L. 



D. The average curvature at local minima 



The next natural question to ask is how the average curvature, K at the minimum points scales with the system size 
for the general system described by eq. (|3^) . This can be evaluated as the conditional average of the local curvature 
at the local minima: 

ir.. u. , ((0.-^._i)e(-^._i)e(0.)) ((^2-^i)9(-0i)9(</)2)) 

A)min = 0* - 0i-l) min = , TTTTXu = ^4) 

(e(-0j_l)e((/)j)) (li) 

where translational invariance is exploited again. The numerator in (^^ can be obtained after performing the same 
basis transformation (Appendix A) that was essential to find (u) . Then after elementary integrations we find 



WvS ^/^S^ ™os(C{(l)/C*(0)) 

Using the explicit results for the slope correlation function for z — 2 and z = 4, and the scaling forms for it for higher 
z given in the previous subsections, one can easily deduce the following. For z < 5 the average curvature at the local 
minimum points on a lattice tends to a constant in the thermodynamic limit. For z = 2 

{KU:n^'-^xf^ + o('-] , (46) 

and for z = 4 



1/ \L 



{KU^c,V^^g+o(^^^ . (47) 
The behavior of this quantity drastically changes for z > 5, where it diverges with the system size as: 

(if),nin ~ (48) 
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III. OTHER LATTICE MODELS AND AN APPLICATION TO PARALLEL COMPUTING 



A. The single-step model 

In the single-step model the height differences (i.e., the local slopes) are restricted to ±1, and the evolution consists 
of particles of height 2 being deposited at the local minima. While the full dynamic behavior of the model belongs 
to the KPZ universality class, in one dimension the steady state is governed by the EW Hamiltonian j2^. Thus, the 
roughness exponent is a = 1/2, and we expect the finite-size effects for (u) to follow eq. (|2^). The advantage of 
this model is that it can be mapped onto a hard-core lattice gas for which the steady-state probability distribution of 
the configurations is known exactly ]25| , p9t . This enables us to find arbitrary moments of the local minimum density 
operator. Since = ±1, it can be simly written as 

7 2^ ^ - 72^(1 -"^-i)'^^ > (49) 



L ^ 2 2 L 

i=l i=l 



where = (1 + (j)i)/2, corresponds to the hard core lattice gas occupation number. The constraint 'Y^^=i4'i — 



translates to X^i^Li = L/2. Note that here {u) — ((1 — ni-i)ni) is proportional to the average current. Knowing 
the exact steady-state probability distribution |23,E3], one can easily find that 



1 L-2 



Thus the exact finite-size effects for the local minimum density: 

1 L 11 



4 L- 1 



(50) 



in qualitative agreement with (| 



B. The Massively Parallel Exponential Update model 



One of the most challenging areas in parallel computing |2^] is the efficient implementation of dynamic Monte- 
Carlo algorithms for discrete-event simulations on massively parallel architectures. As already mentioned in the 
Introduction, it has numerous practical applications ranging from magnetic systems (the discrete events are spin 
fiips) to queueing networks ( the discrete events are job arrivals). A parallel architecture by definition contains 
(usually) a large number of processors, or processing elements (PE-s). During the simulation each processor has 
to tackle only a fraction of the full computing task (e.g., a specific block of spins), and the algorithm has to ensure 
through synchronization that the underlying dynamics is not altered. In a wide range of models the discrete events are 
Poisson arrivals. Since this stochastic process is reproducible (the sum of two Poisson processes is a Poisson process 
again with a new arrival frequency) , the Poisson streams can be simulated simultaneously on each subsystem carried 
by each PE. As a consequence, the simulated time is local and random, incremented by exponentially distributed 
random variables on each PE. However, the algorithm has to ensure that causality accross the boundaries of the 
neighboring blocks is not violated. This requires a comparison between the neighboring simulated times, and waiting, 
if necessary (conservative approach). In the simplest scenario (one site/PE), this means that only those PEs will 
be allowed to attempt the update the state of the underlying site and increment their local time, where the local 
simulated time is a local minimum regarding the full simulated time horizon of the system, {ri}, i = 1,..,L (for 
simplicity we consider a chain-like connectivity among the PE-s but connectivities of higher degree can be treated 
as well). One can in fact think of the time horizon as a fluctuating surface with height variable r^. Other examples 
where the update attempts are independent Poisson arrivals include arriving calls in the wireless cellular network of a 
large metropolitan area j2j, or the spin flip attempts in an Ising ferromagnet. This extremely robust parallel scheme 
was introduced by Lubachevsky, and it is applicable to a wide range of stochastic cellular automata with local 
dynamics where the discrete events are Poisson arrivals. The local random time increments is, in the language of the 
associated surface, equivalent to depositing random amounts of 'material' (with an exponential distribution) at the 
local minima of the surface, see Figure 2. This in fact deflnes a simple surface growth model which we shall refer to as 
'the massively parallel exponential update model' (MPEU). The main concern about a parallel implementation is its 
efficiency. Since in the next time step only a fraction of PE-s will get updated, i.e., those that are in the local minima 
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of the time horizon, while the rest are in idle, the efficiency is nothing but the average number of non-idling PE-s 
divided by the total number of PE-s (L), i.e., the average number of minima per unit length, or the minimum-point 
density, u. The fundamental question of the so called scalability arises: will the efficiency of the algorithm go to zero 
as the number of PE-s is increased (L — > oo) indefinitely, or not? If the efficiency has a non-zero lower bound for 
L ^ OO the algorithm is called scalable, and certainly this is the preferred type of scheme. Can one design in principle 
such efficient algorithms? 




FIG. 2. The MPEU model. The arrows show the local minima where random amounts of material will be deposited in the 
next time step. 

As mentioned in the Introduction, we know of one example that nature provides with an efficient algorithm for a 
very large number of processing elements: the human brain with its 10^"'^ PE-s is the largest parallel computer ever 
built. Although the intuition suggests that indeed there are scalable parallel schemes, it has only been proved recently, 
see for details Ref. §, by using the aforementioned analogy with the simple MPEU surface growth model. While 
the MPEU model exactly mimics the evolution of the simulated time-horizon, it can also be considered as a primitive 
model for ion sputtering of surfaces (etching dynamics): to see this, define a new height variable via hi = — r,;, i.e. flip 
Figure 2 upside down. This means that instead of depositing material we have to take, 'etch', and this has to be done 
at the local maxima of the {hi} surface. In sputtering of surfaces by ion bombardement an incoming ion-projectile 
will most likely 'break off' a piece from the top of a mound instead from a valley, very similar to our 'reversed' 
MPEU model. It was shown that the sputtering process is described by the KPZ equation, |^ , |l^ . This qualitative 
argument is in complete agrrement with the extensive MC simulations and a coarse-grained approximation of Ref. |^ 
that MPEU, similar to the single-step model, it also belongs to the KPZ dynamic universality class; in one dimension 
the macroscopic landscape is governed by the EW Hamiltonian. 

The slope varaibles 0i for MPEU are not independent in the L — > oo limit, but short-ranged. This already guarantees 
that the steady-state behavior is governed by the EW Hamiltonian, and the density of local minima does not vanish 
in the thermodynamic limit. Our results confirm that the finite-size effects for {u) follow eq. (Eq): 



{u)l 



const. 



(52) 



with 



= 0.24641(7), see Fig. 3. 
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FIG. 3. Density of minima vs. 1/L for the MPEU model. 

We conclude that the basic algorithm (one site per PE) is scalable for one-dimensional arrays. The same correspon- 
dence can be applied to model the performance of the algorithm for higher-dimensional logical PE topologies. While 
this will involve the typical difficulties of surface-growth modeling, such as an absence of exact results and very long 
simulation times, it establishes potentially fruitful connections between two traditionally separate research areas. 

C. The larger curvature model 

In this subsection we briefly present a curvature driven SOS surface deposition model known in the literature as 
the larger curvature model, and show a numerical analysis of the density of minima on this model. This model was 
originally introduced by Kim and Das Sarma and Krug ||27| independently, as an atomistic deposition model 
which fully conforms to the behaviour of the continuum fourth order linear MuUins equation (i/ = 0, k > in Eq. |^). 
Note that the discrete analysis we presented in Section II is based on the discretization of the continuum equation 
using the simplest forward Euler differencing scheme. The larger curvature model, however, is a growth model where 
the freshly deposited particles diffuse on the surface according to the rules of the model until they are embedded. 
Since in all the quantities studied so far, the correspondence (on the level of scaling) between the larger curvature 
model and the Langevin equation is very good, we would expect that the dynamic scaling properties of the density 
of minima for both the model and the equation to be identical. 

The large curvature model has rather simple rules: a freshly deposited atom (let us say at site i) will be incorporated 
at the nearest neighbor site which has largest curvature (i.e., Ki — /i^+i + /i^-i — 2ft,,; is maximum). If there are more 
neighbors with the same maximum curvature, then one is chosen randomly. If the original site (i) is among those 
with maximum curvature, then the atom is incorporated at i. 

Figure 4 shows the scahng of the density of minima < u >l in the steady state, vs. 1/\/L. According to Eq. (|3^), 
for the fourth order equation on the lattice, the behavior of the density of minima in the steady state scales with 
system size as I/VL- And indeed. Figure 4 shows the same behaviour for the larger curvature model, as expected. 
Note that this behavior sets in at rather small system sizes already, at about L — 100, meaning that the finite system 
size effects are rather small for the larger curvature model. This is a very fortunate property since increasing the 
system size means decreasing the density of minima, therefore relative statistical errors will increase. 
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FIG. 4. Density of minima in the steady state for the larger curvature model. 

This can only be improved by better statistics, i.e., with averages over larger number of runs. This becomes however 
quickly a daunting task, since the cross-over time toward the steady state scales with system size as L'^. As we shall 
see in Section V.A, a matematically rigorous approach to the continuum equation yields the same behaviour. 
Since the density of minima does decay to zero, an algorithm corresponding to the larger curvature model (or the 
MuUins equation) would not be asymptotically scalable. 

Finally, we would like to make a brief note about the observed morphologies in the steady state for the MuUins 
equation , or the related models. It has been shown previously that in the steady state the morphology tipically 
shows a single large mound (or macroscopic groove). At first sight this may appear as a surprise, since we have 
shown that the number of minima (or maxima) diverges as ^/L (the density vanishes as 1/\/L- There is however 
no contradiction, because that refers to a a mound that expands throughout the system, i.e. it is a long wavelenght 
structure, whereas the number of minima measures all the minima, and thus it is a short wavelenght characteristic. 
In the steady state we indeed have a single large, macroscopic groove, however, there are numerous small dips and 
humps generated by the constant coupling to the noise. 

IV. EXTREMAL-POINT DENSITIES ON THE CONTINUUM 

Let us consider a continuous and at least two times differentiable function / : [0, L] R. We are interested in 
counting the total number of extrema of / in the [0, L] interval. The topology of continuous curves in one dimension 
allows for three possibilities on the nature of a point Xi for which f'\xi = 0. Namely, Xi is a local minimum if f"\xi > 0, 
a local maximum if f"\xi < and it is degenerate if f"\xi — 0. We call the point Xi a degenerate flat of order k, if 
/^"'■'Ui = for j = 1,2, ..,k and f'^'^^'^^xi 7^ 0, fc > 2, assuming that the higher order derivatives /'^■'^ implied exist. 
The counter-like quantity 

L 

c{L,[f])^^Jdx\f"\d{n (53) 


where S is the Dirac-delta, gives the number of extremum points per unit length in the interval [0,1/], which in the 
limit of L ^ is the extremum point density of / in the origin. For our purposes L will always be a finite number, 
however, for the sake of briefness we shall refer to c simply as the density of extrema. Note that counting the extrema 
of a function / is equivalent to counting the zeros of its derivate /'. The divergence of c for finite L implies either the 
existence of completely flat regions (infinitely degenerate), or an "infinitely wrinkled" region, such as for the truncated 
Weierstrass function shown in Fig. 2 ( in this latter case the divergence is understood by taking the limit M 00). 
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As already explained in the Introduction this infinitely wrinkled region does not necessarily imply that the curve is 
fractal, but if the curve is fractal, then regions of infinite wrinkledness must exist. The divergence or non-divergence 
of c can be used as an indicator of the existence of such regions (completely flat or infinitely wrinkled). 

One can make the following precise statement related to the counter c: if Xi is an extremum point of / of at most 
finite degeneracy fc, and if there exist a small enough e, such that / is analytic in the neighborhood [xi — e, Xi + e], 
and there are no other cxtrcma in this neighborhood, then 



/(xO= J dx\f"\Sif') = l, 0<e«l 



(54) 



In the following we give a proof to this statement. 
Using Taylor-series expansions around Xi, one writes: 



nx)^^{x-x.)'^ + 



ki 



(fcTT) 



■{x - x,)''+' + ... 



5{g{x))=Y^\g\x,r^ 5{a 



Xi are simple zeros of g. 



(55) 
(56) 



where we introduced the shorthand notation aj = f'^-'^'^^\xi- For the non-degenerate case of fc = 1, ( p4[ ) follows from 
a classical property of the delta function, namely: 



(57) 



Let us now assume that Xi is degenerate of order k (k > 2). Using the expansions (|55|), (pq), the variable change 
u = X — Xi, and the well-known property 5(ax) = \a\~^6{x), we obtain: 



I{xi) ^ k du \u\ 



^ (k-i 



jk - 1)! ak+j 



-(fc-l+i)! ak 



1 + E 



A:! 



^ {k + jy. ak 



u-' 



(58) 



Next we split the integral ( |58| ) in two: ^ ... — ... + JJ, make the variable change u —u in the first one, and 
then z = m'^ in both integrals. The final expression can then be written in the form: 



lix,) = / dz \A{z)\ 5{zB{z)) , 



(59) 



where 



{k - 1)! ak+j 



, and B{z) = 1 + ^ 



k\ ak+ 



j-{{k + j)\ ak 



j^^{k-l+j)\ ak 
We have A(0) = B(0) = 1, and 



i zHzli-i 



.B(z)r = l + ^ 



^(fc + j)! ak 



(60) 



(61) 



(Take the derivatives separately to the right and to the left of z = 0). Thus, since z = is a simple zero of zB{z) 
property (^7|) can be applied for sufficiently small e: 



I{x,) ^ |A(0)| = 1 



(62) 



proving our assertion. Note that because of (54), c counts all the non-degenerate and the finitely degenerate points 
as well, giving the equal weight of unity to each. Can we count the non-degenerate extrema separately? The answer 
is affirmative, if one considers instead of (p3) the following quantity: 



c,(L,[/])^i / dx\rr^5{r), q>o 



(63) 
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Performing the same steps as above we obtain for a degenerate point: 

Xi+e 



(x,) EE I \fr+' Sif) = J dz \z\<'-^) \A{z)r^ 5{zB{z)) . (64) 



Since fc > 2, g (l - i) > > 0, i.e., 

Iq{xi) ^ , for Xi degenerate. (65) 

This means, that g > ehminates the degenerate points from the count. To non-degenerate points (fc = 1) (|63| ) gives 
the weight of 

q 

Iq{xi) — \ai\'^ — f"\xi , for Xi non-degenerate. (66) 

In other words, 

Cq{L, [/]) = ^ ^ l^i^i)]"^ J q > 0, Xi non-degenerate extrema of / (67) 

i 

where K{x) = f" is the curvature of / at x. The hmit g ^ 0+ in ( |67[ ) gives the extremum point density c(L, [/]) of 
/ of non-degenerate extrema: 

L 

c(L, [/]) = hm Cq{L, [/]) = lim j f dx |/'T+^ S{f) (68) 

g— »0+ g— ♦0+ Li J 



It is important to note, that taking the g — > 0^ limit in (|67| ) is not equivalent to taking q = in (|63|), i.e., the limit 
and the integral on the rhs of (^8|) are not interchangeable! The difference is the set of degenerate points! 

Until now, we did not make any distiction between maxima and minima. In a natural way, we expect that the 
quantity: 

L 

u{L, [/]) dx f" Sif) eif") (69) 



where 9{x) is the Heaviside step-function, will give the density of minima (due to the step function, here we can 
drop the absolute values). However, performing a similar derivation as above, one concludes that (69) is a little bit 
ill-defined, in the sense that the weight given to degenerate points depends on the definition of the step-function in 
the origin (however, u{L, [/]) is bounded). Introducing a, q — regulator as above, the weight of degenerate points is 
pulled down to zero: 

L 

Uq{L, [/]) ^\jdx eif") , q>0. (70) 



and 

Uq{L, [/]) = — [isr(xi)]^ , q > 0, Xi non-degenerate minima of / (71) 

i 

Note that in the equation above the absolute values are not needed, since we are summing over the curvatures of all 
local minima. The density u{L, [/]) of non-degenerate minima of / is obtained after taking the limit g — > 0+: 

u{L,[f])= Ih-n Uq{L,[f]) (72) 

and the limit is not interchangeable with the integral in (|7^) . To obtain densities for maxima, one only has to replace 
the argument /" of the Heaviside function with — /". 
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A. Stochastic extremal-point densities 



We are interested to explore the previously introduced quantities for a stochastic function, subject to time evolution, 
h{x,t). This function may be for example the solution to a Langevin equation. We define the two basic quantities in 
the same way as before, except that now one performs a stochastic average over the noise, as well: 



dx 



9+1 



^1? 

OX 



and Uq{L,t) 



dx^ 



9+1 



ox J \ Ox'^ J 



(73) 



For systems preserving translational invariance, the stochastic average of the integrand becomes a;- independent, and 
the integrals can be dropped: 



C,(i,t) 



C/,(L,t) 



dx"^ 
dx"^ 



9+1 



9+1 



dx ^ 



dx'^ 



(74) 
(75) 



According to ( pT\ j and (71), Cq{L,t) and Uq{L,t) can be thought of as time dependent "partition functions" for 
the non-degenerate extremal-point densities of the underlying stochastic process, with q playing the role of "inverse 
temperature" : 



Cq{L, t) — Iy ^ \K{xi)\'^ \ , q> Xi non-degenerate extrema 



(76) 



Uq{L,t) — / — [K{xi)]''\ , q > 0, Xi non-degenerate minima 



(77) 



It is important to mention that in the above equations the average < .. > and the summation arc not interchangeable: 
particular realizations of h have particular sets of minima. 

Two values for q are of special interest: when g — > 0"*" and q = I. In the first case we obtain the stochastic average 
of the density of non-degenerate extrema and minima: 



C(L,i)= lim Cq{L,t) , and U{L,t) = lim Uq{L,t) , 

♦0+ q^0+ 



and in the second case we obtain the stochastic average of the mean curvature at extrema and minima: 

and Kmm y 



C{L,t) 



I) — -=- 



U{L,t) 



(78) 



(79) 



(we need to normalize with the number of extrema/minima per unit length to get the curvature per ex- 
tremum/minimum). 

In the following we explore the quantities (|7^)-(|79|) for a large class of linear Langevin equations. To simplify the 
calculations, we will assume that g is a positive integer. Then we will attempt analytic continuation on the final result 
as a function of q. In the calculations we will make extensive use of the standard integral representations of the delta 
and step functions: 



dz 



00 



"(!/) = li-n / ^ ' 



n=0_ 
izy 



dz {izY' 



In n\ 



n—O 



1 



27r e + iz n\ 



(80) 



(81) 



If g is a positive integer, we may drop the absolute value signs in (|75|). In ( |74| ) we can only do that for odd q. 
The absolute values make the calculation of stochastic averages very difficult. We can get around this problem by 
employing the following identity: 



|yr = j/"{(-i)" + %) [i-(-i)"]} 



(82) 
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This brings to 



C,{L,t) = [1 - (-!)«] Ug{L,t) + {-iy+'B,{L,t) (83) 

where 

9+1 I /a2u\ 9+1 



Obviously for q odd integer, Bq = Cq. For g even, Bq is an interesting quantity by itself. In this case the weight of an 
extremum Xi is sgn{K{xi))\K{xi)\'^ . If the analytic continuation can be performed, then the q —> Q'^ limit will tell us 
if there are more non-degenerate maxima than minima (or otherwise) in average. Using the integral representations 
(|^) and (|§): 

°° dz (izY I fd^hV^^ /dh' ' 
2n nl \\dx^) \dx ^ 

f7,(L,t) = ^hm ^ ^ y ——^y ——P-^ n^j / 

ni=0n2=0_^ _„ \ ^ ' ^ ' I 



n=0 



V. EXTREMAL-POINT DENSITIES OF LINEAR STOCHASTIC EVOLUTION EQUATIONS 

Next we calculate the densities ( |85| ) , ( p6| ) for the following type of linear stochastic equations: 

^^-iy{-W^y^^h + T]{x,t) , iy,D,z>0, xe[0,L] (87) 

with initial condition h{x,0) ~ 0, for all a; G [0, i]. ?7 is a white noise term drawn from a Gaussian distribution with 
zero mean {r]{x, t)) = 0, and covariance : 

{rj{x, t)ri{x',t')) = 2D5{x - x') 5{t - t') , (88) 

We also performed our calculations with other noise types, such as volume conserving and long-range correlated, 
however the details are to lengthy to be included in the present paper, it will be the subject of a future publication. 
As boundary condition we choose periodic boundaries: 

h{x + nL, t) = h{x, t) , ri{x + nL, t) = ri{x, t) , for all n G Z (89) 



The general solution to (g7|) is obtained simply with the help of Fourier series |18|. The Fourier series and its 
coefficients for a function / defined on [0, L] is 

L 



where k = ^n, n = .., —2, —1, 0, 1, 2, ... The Fourier coefficients of the general solution to (87) are: 

t 

h{k, t) = j dt' e-''l'=l'(*-*')?7(fc, t') (91) 



The correlations of the noise in momentum space are: 

2D 

{fj{k, t)^{k', t')) = — Sk^-k' Sit - t') . (92) 
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Due to the Gaussian character of the noise, the two-point correlation of the solution (|91|) is also delta-correlated 
and it completely characterizes the statistical properties of the stochastic dynamics (87). It is given by: 



h{k,t)h{k\t')) ^ S{k,t)5k^^u' 



(93) 



where S(k,t) is the structure factor given by: 



5(.,t) = -^[l-e-l'=l-] . (94) 

Equation (|^ has been analyzed in great detail by a number of authors, see Ref. 11 for a review. It was shown that 
there exist un upper critical dimension dc = z for the noisy case of Eq. (^) which separates the rough regime with 
d < z from the non- roughening regime d > z. In one dimension, the rough regime corresponds to the condition z > 1, 
which we shall assume from now on, since this is where the interesting physics lies. 

Next, we evaluate the quantities ([7^-([79|) via directly calculating the expressions in ( |85| ) and (|6|). This amounts 
to computing averages of type: 



Expressing h with its Fourier series according to (|90|), we write: 

-j =i*^^...^fci...A:M /i(fci,i).../i(fcAf,i) e^('=i+-+'=")" (96) 

ki kM 

= i-^f E - E k[\..k'^^ h{k[,t)..:h{k'N,t) e^«+-+'=«)- (97) 



,2, N N 



dx 

which then is inserted in (^5|). Thus in Fourier space one needs to calculate averages of type 
{h{ki,t)...h{kM,t)h{k'i,t)...h{k'j^A)). According to ( p3| ) h is anti-delta-correlated, therefore these averages can be 
performed in the standard way |pO| which is by taking all the possible pairings of indices and employing (|9^). In our 
case there are three types of pairings: {kj, ki}, {kj, fc;}, and {k'j, fc;}. Let us pick a 'mixed' pair {kj, k'l} containing a 
primed and a non-primed index. The corresponding contribution in the Qn,m will be: 

E ''A'sm, t)e^^'^+'>6,^^_,, (98) 

Since the structure factor S{k,t) is an even function in k, ( p8| ) becomes X^fe^ f^j^S!{kj,t) = 0, because the summand 
is an odd function of kj and the summation is symmetric around zero. Thus, it is enough to consider non-mixed 
index-pairs, only. This means, that Qn,m decouples into: 



;n,m = 




(99) 



The averages are calculated easily, and we find: 

I fdh\^^\ { (^^-1)!! [(t>2{L,t)]^''^ , for M even , 



dx 



(100) 



, for M odd 



and 

^\ { [N [(t)i{L,t)]'^'^ , for even 



(101) 



, for M odd 
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where 



(102) 



Employing ( |100| ), and (101) in (p5[), it follows that if q is an even integer, q — 2s, s = 1, 2, ..: 

52.(0=0, s = l,2,... 
whereas for q odd integer, q = 2s — 1, s = 1, 2, ..: 



B2s-l{t) — C2s-l(i) 



2'*-2 



r u 



1\ [^i{L,t)]' 



2/ y^^CM) 



, s = l,2,... 



(103) 



(104) 



where we used the identity 2^(2^— l)!!/(2p)! — 1/p!, and performed the Gaussian integral. 

The calculation of Uq is a bit trickier. The sum over ni in ( |8^ ) is easy and leads to the Gaussian e~'^2(i,t)zi/2 
However, the sum over n2 is more involved. Let us make the temporary notation for the sum over n2- 



= E (2r - 1)!! [04(i, t)]'' , n2 + g + 1 = 2r 



n2— 



We have to distinguish two cases according to the parity of q: 
1) q is odd, g = 2s — 1, s = 1, 2, ... In this case Rq becomes 



R2S-1 = ^ 



[2{r-s)]\ r\ 

The Hermite polynomials are defined via the Rodrigues formula as: 



Qj.2s 



x=l 



77„(x)^(-l)"e-~33^(e--') 



Using this, we can express i?2s-i with the help of Hermite polynomials: 

'1 



R2S-1 = i-iy 



2s 



2) q is even, g = 2s, s = 1, 2, ... The calculations are analogous to the odd case: 



{iZ2) 



^^_^^[2{r-s)-l]\ r! 
or via Hermite polynomials: 



(iZ2)~ 



Q^2s+l 



-4>4,{L.t)zlx^ /2 



R2s = i{-iy 



H2S+1 [ \l-ML,t)z2\ e-*-(^'*)^^^-^/2 



(105) 



(106) 



(107) 



(108) 



(109) 



(110) 



In order to obtain Uq we have to do the integral over Z2 in (B6|). This can be obtained after using the formula: 



dx {x ± icY Hn{x) e 



r (2— iil 

M 2 ; g±^(.+n) c^0+. 



r(-i.) 



(Ill) 



Finally, the densities for the minima read as: 



^.(0^^r(s + i) 



U2s-l{t) 



r s 



V 2; V0;(z;i) 



(112) 
(113) 
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Formulas (103), (104), (112), and (113) combined with (83) can be condensed very simply, and we obtain the general 
result as: 



U,iL,t) 



TT 

2Uq{L,t) 



(I-') 



[<j^A{L,t)] — 



(114) 
(115) 



Equa tion s (114), (115) with together with (113) fully solve the problem for the density of non-degenerate extrema. 
Eq. (|115| ) is an expected result in one dimension, because Eq. ( pTj ) preserves the up-down symmetry. The density of 
non-degenerate minima is: 



and the stochastic average of the mean curvature at a minimum point is: 



U{L,t) 



(116) 



(117) 



i.e., the average curvature at a minimum is proportional to the square root of the fourth moment of the structure 
factor. In the following section we exploit the physical information behind the above expressions for the stochastic 
process (^7|). At some parameter values a few, or all the quantities above may diverge. In this case we introduce a 
microscopic lattice cut-off < a <C 1, and analyze the limit a — > 0"*" in the final formulas. This in fact corresponds 
to placing the whole problem on a lattice with lattice constant a. It has been shown in Ref. ||l^ that for the class 
of equations (87) there are three important length-scales that govern the statistical behavior of the interface h: the 
lattice constant a, the system size L, and the dynamical correlation length ^ defined by: 



According to (|102) and (M) the function 0„i(L,t) becomes: 



(pm{L,t) = 



27rn 

IT 



/ 2-Kn 

= L 



m = 2,4 



(118) 



(119) 



The n — Q term can be dropped from the sum above, because it is zero even for m < z (expand the exponential 
and then take n — 0). However, the whole sum may diverge depending on m and z. In order to handle all the cases, 
including the divergent ones we introduce a microscopic lattice cut-off a, < a ^ 1, and then analyze the limit 
a — !■ 0"*" in the fin al ex pressions. This is in fact equivalent to putting the whole problem on a lattice of lattice spacing 
a. Appropiately, (119) becomes: 



27m ^ 
L ) 



m = 2,4 



(120) 



A. Steady-state regime. 



Putting ^ = CX3 in (120) takes a simpler form: 
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7n — Z 2a 



m = 2,4 



(121) 



As a ^ 0+, 0m becomes proportional to C(z — m). For z — m > 1 (j)m is convergent, otherwise it is divergent. In the 
divergent case we quote the following results: 



AT 



= \nN + C + Oil/N) , if s = -1 



(122) 
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and 



N 



S+1 



1 + Oil/N) 



if s > -1 



(123) 



which we will use to derive the leading behaviour of the extremal point densities when L/a oo. From equations 
(Iml), (PI), (pl6|) and (pTzl) follows: 



2D 





"L/2a 


9+1 
2 


"L/2a 


(27r)*(5-^)L-i-«(5-) 


n=l 







(124) 
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and 



K{L, oo) 
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L/2a 



(126) 



The convergency (divergency) properties of the sums in Eqs. ( 124 - 126 ) for a — > 0^ generate two critical values for z, 
namely z — i and z = 5. In the three regions separated by these values we obtain qualitatively different behaviors for 
the extremal-point densities. 

i) z > 5. All quantities are convergent as a 0+. We have: 



U,iL, oo) = r f I + l) f (2.)f i-i+f (-5) 

' ' V2 y \nuj ^ ' [C(^-2)]i 



z > 5 



(127) 



K{L, oo) = (27r)' 





-4) 




-2) ' 



C{z-4)L' 
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(129) 



Eq. ( |128| ) shows that there are a finite number of minima {\/C{^ ~ 4)/C(z — 2)) in the steady state, independently of 
the system size L. {U{L, oo) is the number of minima per unit length, and LU (L, oo) is the number of minima on the 
substrate of size L). The mean curvature K{L, oo) diverges with system size as L^^"^^/^. This is consistent with the 
fact that the system size grows as L, the width grows as L^^~^^^'^, i.e., faster than L, and thus the peaks and minima 
should become sleeker and sharper as L — > oo, expecting diverging curvatures in minima and maxima. However, this 
is not always true, since the sleekness of the humps and mounds does not necessarily imply large curvatures in minima 
and maxima if the shape of the humps also changes as L changes, i.e., there is lack of self- affinity. The existence of 
z = 5 as a critical value is a no n-trivial results coming from the presented analysis. 

a) z = 5. According to (122), <j)4{L, oo) diverges logarithmically as a ^ 0+. One obtains: 



c/,(oo) r (I + 1) 
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(130) 



U{L, oo) = ^Jln^+C 
^ ' ^C(3) i V 2a 



(131) 
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K{L, oo) = J-^ [\n^+C 



2v 



2a 



(132) 



Eq. (131) shows that the although the density of minima vanishes, the number of minima is no longer a constant but 
diverges logarithmically with sy stem size L. The mean curvature still diverges, but logarithmically, when compared 

to the power law divergence of (12£)- 

For the mean curvature K(oo) in ( |126| ) z = 5 is the only critical value, since it only depends on (l>4. For z < 5, 
using Eq. (123) we arrive to the result that the mean curvature in a minimum point approaches to an L-independent 
constant for L/a oo with corrections on the order of a/L: 



D 



2v{b - z) 



z < 5 



(133) 



We arrived to the same conclusion in Section II. D when we studied the steady state of the discretized version of the 
continuum equation. Coincidentally, for z = 4 the two constant values from (133) and ( ^7| ) are identical (a = 1 by 
definition in (flTj). 

Hi) 3 < z < 5. In this case 4'a{L^ oo) oo and 4>2{L, oo) < oo as a — > 0+, and: 
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.'^-'-^(l-)(S)'(^ 



|(5-z) 



(5-z)^v'C(^ 



and 



U{L, oo) 



v/(5 - z)C(z - 2) 



(134) 



(135) 



and the mean curvature is just given by (|l33|) . 

Comparing Eqs. (127), (130), and ( [134| ) we can make an interesting observation: while for z > 5 the dependence 
on the system size L is coupled to the 'inverse temperature' q, for 3 < z < 5 the dependence on L decouples from g, 
i.e., it becomes independent of the inverse temperature! Eq. (135;) shows that the density of minima vanishes with 
system size as a power law with an exponent (z — 3)/2 but the number of minima of the substrate diverges as a power 
law with an exponent of (5 — z)/2. 

iv) z = 3. In this case (/>4(i, oo) — > oo and 02(i, oo) — > oo logarithmically as a — s- 0^. One obtains: 
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(136) 



and 
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2V2a 



(137) 



with a logarithmically vanishing density of minima, and the dependence on the system size in ( |136| ) is not coupled 
to q. 



v) 1 < z < 3. Now both (/)4 and (j>2 diverge as a ^ 0+. Employing (123), yields: 

'2D\^ /7r\f(5-2) 
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Note, that in leading order, both Uq{L, oo) and the density of minima U{L, oo) become system size independent! 
The system size dependence comes in as corrections on the order of a/L and higher. The fact that the efficiency of 
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the massively parallel algorithm presented in Section III.B is not vanishing is due precisely to the above phenomenon: 
the fluctuations of the time horizon in the steady state belong to the z = 2 class (Edwards- Wilkinson universality), 
and according to the results under iv), the density of minima (or the efficiency of the parallel algorithm) converges 
to a non-zero constant, as L —^ oo, ensuring the scalability of the algorithm. An algorithm that would map into a 
z > 3 class would have a va nishing efficiency with increasing the number of processing elements. In particular, for 
z = 2, one obtains from ([l39| ) U{L, oo) ~ (a2\/3)"i = 0.2886. ../a. Note that the utilization we obtained is somewhat 
different from the discrete case which was 0.25. This is due to the fact that this number is non-universal and it may 
show differences depending on the discretization scheme used, however it cannot be zero. 

Another important conclusion can be drawn from the final results enlisted above: at and below z = 5, all the 
quantities diverge when a — > 0"*", and keep L fixed. This means that the higher the resolution the more details we find 
in the morphology, just as for an infinitely wrinkled, or a fractal-like surface. We call this transition accross z = 5 a 
'wrinkle' transition. As shown in the Introduction, wrinkledness can assume two phases depending on whether the 
curve is a fractal or not and the transition between these two pases may be conceived as a phase transition. However, 
one may be able to scale the system size L with a such that the quantities calculated will not diverge in this limit. 
This is possible only in the regime 3 < z < 5, when we impose: 



-3Q(g+l)(5-2:) ^ 



= const. 



This shows that the rescaling cannot be done for all inverse temperatures q at the same time, 
density of minima and z = 4, La = const. 



(140) 

In particular, for the 



B. Scaling regime 



In order to obtain the temporal behavior of the extremal-point densities we will use the Poisson summation formula 
(B4) from Appendix B on (|12C|). After simple changes of variables in the integrals this leads to: 
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x^~^ cos ( jnx 1(1- e-^' 



(141) 



This expression shows, that the scaling properties of the dynamics are determined by the dimensionless ratios i/^ 
and ^/a. The scaling regime is defined by a ^ ^ <C i. 

As we have seen in the previous section, <j)m is convergent for z > m + 1 but diverges when z<m-|-l, asa^O^. 



In the convergent case, the lattice spacing a can be taken as zero, and thus the first term on t he r hs. of (141) vanishes 
and the time dependence of the infinite system-size piece of (j)m (the first integral term in ( 141 )) assumes the clean 
po wer-l aw behaviour of ^l^-™-!)/^ (with a positive exponent). In the divergent case, however, the non-integral term 
of (141) does not vanish, and the time-dependence will not be a clean power-law. Even the integral terms will present 
corrections to the power- law (which has now a negative exponent), since the limits for integration contain 

^. The first intergal on the rhs of (141) for z m + 1 can be calculated exactly: 
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(142) 



where T{a,x) is the incomplete Gamma function. In our case (ir^/a)^ is a large number, and therefore we can use 
the asymptotic representation of T{a,x) for large x, see pp. 951, equation 8.357. According to this, for large 
x, T(a,x) ~ a;° ~^e~ ^, i.e., it can become arbitrarily small, with an exponential decay. Thi s ter m can the refore be 

we will 



neglected from (142), compared to the othe r two terms, even in the divergent case. Inserting ( |l42| ) into ( |l4l[ ) 
see that also the non- integral piece of ( 141 ) can be neglected compared to the term generated by the first on the rhs 
of ( |142| ), since in the scaling regime a <C ^ ^ i, and thus the ratio a/L can be neglected compared to (m — z + 
(This is needed only in the divergent regime, z < m + 1.) Thus, one obtains: 
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(143) 
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where 
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oo t. 

-jjp- ^ / a;™"^ cos (Anx) (^1 - e""^"^ , z 7^ m + 1 



(144) 



The oscillating terms condensed in will give the finite-size corrections, as long as L/^ 3> 1. 
The 2; = m + 1 case (divergent) can also be calculated, however, instead of (142) now we have: 
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i_,-n=i„i^iU2-iEir-ri^ 



, z — m + 1 



(145) 



where Ei(x) is the exponential integral function. According to the large-x expansion of the exponential integral 
function, see |Q, pp. 935, equation 8.215, Ei(— a;) ~ — a;~^e~^, it is vanishing exponentially fast, thus it can be 
neglected in the expression of (f>m in the scaling limit: 
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where 
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Fm{\, p) = 2 ^ j — cos {Xnx) [l 



, Z = TO + 1 



(147) 



Thus in the scaling limit, the temporal behavior of becomes a logarithmic time dependence plus a constant, as 
long as i/^ » 1. 

Observe that for z < m + 1 t he f irst term on the rhs. of ( 143 ), reprod uces exa ctly the diverging term (as a — s- 0^) 
of the steady-state expression (121) which can be seen after employing ( |122| ) in (121). This means that for ^ — > 00, 
EmjL /£,,£,/ a) dive rges slower than ^^j^jg jg j^q.^^ ^]-^g saturation occurs). Similarly, for z = m-|- 1 the first term 

on th e rhs of ( |146| ) (after replacing ^ with L ) r epro duce s the diverging term (as a — > 0+) of the steady-state expression 



(121) which can be seen af ter e mploying (123) in ( |l2l|) . This means, that in the saturation (or steady-state) regime 
the remaining terms from ( |146| ) must behave as const. + 0{a/ L) + In (L/^), as ^ — > 00 while keeping L and a fixed. 

Just as in the case of steady-state one has to distinguish 5 cases depending on the values of z, with respect to the 
critical values 3 and 5. For the sake of simplicity of writing, we will omit the arguments of Em{X, p) and Fm{X,p)- 

i) z > 5. We have: 
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From Eqs. (114), mm and (117), it follows 
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(149) 



and 



i^(^,o = \/^^K(0]-yr3i^ (151) 

and therefore the time-behaviour is a clean power- law: Uq{L,t) decays as ~ i-[2+9(2-5)]/2^^ U{L,t) - t-'^/\ and 
K{L,t) diverges as - t(^-5)/2^, for L/S, > 1. 
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a) z = 5. In this case (j)/^ takes the form (146) but (j)2 is still given by (|143|). The quantities of interest become: 
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One can observe that the leading temporal behaviour has logarithmic component due to the borderline situation: 
Uq{L,t) decays as ~ t-i/5(int)(9+i)/2^ UiL,t) ~ t-^^^{lnty/^ , and K{L,t) diverges as ~ (Int)!/^. 
in) 3 < z < 5. 
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Density of minima for the larger curvature model as a function of time (the nr of deposited layers) 
100 (diamonds) and L = 120 (crosses). The straight line corresponds to the behavior t~^^^. 



for two system 
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An important conclusion that can be drawn from these expressions is that below z = 5, the leading time-dependence 
of the partition function Uq{L,t) becomes independent oi the inverse temperature q and it presents a clean power-law 
decay ~ t^i^-'^)/'^^ which is the same also for U{L,t). In particular, for 2 = 4 this means a i^^/^ deca y wh ich is 
very well verified by the larger curvature model from Section III.C, see Figure 5. Also notice from Eq. ( |156| ) that 
the leading term is system size independent. And indeed, this property is also in a very good agreement with the 
numerics on the larger curvature model from Figure 5, where the two data sets for L = 100 and L = 120 practically 
coincide. 

Since the mean curvature depends on 04, only, for all cases below z = 5 the dependence is given by the same 
formula (157) (just need to replace the corresponding value for z). 

iv) z = 3. This is another borderline situation, the corresponding expressions are found easily: 
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and the leading time dependences are: Uq{L,t) ~ (Ini) U{L,t) ~ (hit) 
v)l< z <3. 
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In this case the partition function and the density of minima all converge to a constant which in leading order is 
independent of the system size. The density of minima was shown in Section II to have this property in the steady- 
state. Here we see not only that but also the fact that all g-moments show the same behavior, and even more, the 
time behavior before reaching the steady-state constant is not a clean power-law, but rather a decaying correction in 
the approach to this constant. The leading term in the temporal correction is of ~ and the next-to-leading 

has ~ t-(5-^)/^ 



VI. CONCLUSIONS AND OUTLOOK 



In summary, based on the analytical results presented, a short wavelength based analysis of interface fluctuations can 
provide us with novel type of information and give an alternative description of surface morphologies. This analysis 
gives a more detailed characterization and can be used to distinguish interfaces that are 'fuzzy' from those that 
locally appear to be smooth, and the central quantities, the extremal-point densities are numerically and analytically 
accessible. The partition function-like formalism enables us to access a wide range of g-momenta of the local curvatures 
distribution. In the case of the stochastic evolution equations studied we could exactly relate these q-momenta to the 
structure function of the process via the simple quantities 4>2 and (^4. The wide spectrum of results accessed through 
this technique shows the richness of short wavelength physics. This physics is there, and the long wavelength approach 
just simply cannot reproduce it, but instead may suggest an oversimplified picture of the reality. For example, the 
MPEU model has been shown to belong in the steady state to the EW universality class, however, it cannot be 
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described exactly by the EW equation in all respects, not even in the steady-state! For example, the utihzation (or 
density of minima) of the MPEU model is 0.24641 which for the EW model on a lattice is 0.25. Also, if one just 
simply looks at the steady-state configuration, one observes high skewness for the MPEU model , whereas the EW 
is completely up-down symmetric. This can also be shown by comparing the calculated two-slope correlators. For a 
number of models that belong to the KPZ equation universality class, this broken-symmetry property vs. the EW 
case has been extensively investigated by Neergard and den Nijs |^^. The difference on the short wavelength scale 
between two models that otherwise belong to the same universality class lies in the existence of irrelevant operators (in 
the RG sense). Although these operators do not change universal properties, the quantities associated with them may 
be of very practical interest. The parallel computing example shows that the fundamental question of algorithmic 
scalability is answered based on the fact that the simulated time horizon in the steady state belongs to the EW 
universality class, thus it has a finite density of local minima. The actual value of the density of local minima in 
the thermodynamic limit, however, strongly depends on the details of the microscopies, which in principle can be 
described in terms of irrelevant operators |Q . 

The extremal-point densities introduced in the present paper may actually have a broader application than stochastic 
surface fluctuations. The main geometrical characterization of fractal curves is based on the construction of their 
Haussdorff-Besikovich dimension, or the 'box-counting' dimension: one covers the set with small boxes of linear size e 
and then track the divergence of the number of boxes needed to cover in a minimal way the whole set as e is lowered 
to zero. For example, a smooth line in the plane has a dimension of unity, but the Weierstrass curve of (^ has a 
dimension of In 6/ In a (for b > a). The actual length of a fractal curve whose dimension is larger than unity will 
diverge when e 0+. The total length at a given resolution e is a global property of the fractal, it does not tell us 
about the way 'it curves'. The novel measure we propose in (^) is meant to characterize the distribution of a local 
property of the curve, its bending which in turn is a measure of the curve's wrinkledness. For simplicity we formulated 
it for functions, i.e., for curves which are single- valued in a certain direction. This can be remedied and generalized by 
introducing a parametrization 7 € [0, 1] of the curve, and then plotting the local curvature vs. this parameter K{'-^). 
The plot will be a single valued function on which now (||) is easily defined. 

Other desirable extensions of the present technique are: 1) to include a statistical description of the degeneracies 
of higher order, and 2) to repeat the analysis for higher (such as d = 2) substrate dimensions. The latter is promising 
an even richer spectrum of novelties, since in higher dimensions there is a plethora of singular points (V/ = 0) 
which are classified by the eigenvalues of the Hessian matrix of the function in the singular point. Deciphering the 
statistical behaviour of these various singularities for randomly evolving surfaces is an interesting challenge. The 
studies performed by Kondev and Henley p4| on the distribution of contours on random Gaussian surfaces should 
come to a good aid in achieving this goal. In particular we may find the method developed here useful in studying 
the spin-glass ground state, and the spin-glass transition problem. And at last but not the least, we invite the reader 
to consider instead of the Langevin equations studied here, noisy wave equations, with a second derivative of the time 
component, or other stochastic evolution equations. 
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APPENDIX A: {Q(-X-i)Q{X2)) FOR GENERAL COUPLED GAUSSIAN VARIABLES 

The expression we derive in this appendix, despite its simplicity, is probably the most important one concerning the 
extremal-point densities of one-dimensional Gaussian interfaces on a lattice. If the correlation matrix for two possibly 
coupled Gaussian variables is given by 

{x\) = {x\)=d>G (Al) 

{X1X2) = c 

then the distribution follows as 
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where V = (P ~ > 0. We aim to find the average of the stochastic variable u = (d{—xi)S{x2)' 

/ dcick2e{-xi)Q{x2)P{xuX2) (A3) 

-oo ^ — oo 

which is simply the total weight of the density P{xi,X2) in the < 0, 0:2 > quadrant. If c = 0, the density is 
isotropic, and (u) = 1/4. In the general case it is convenient to find a new set of basis vectors, where the probability 
density is isotropic (of course the shape of the original quadrant will transform accordingly). Introducing the following 
linear transformation 

2^2 = \/ -7 



2 V Vd + c y/d — c 
and exploiting that Q{\x) — Q{x) for A > we have 
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Now the probability density for the new variables, yi,?/2, is isotropic, and (u) — 9/{2'k), where 9 is the angle enclosed 
by the following two unit vectors: 
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From their dot product one obtains 



and, thus, for (u) 



cos{e) = = ^. (A7) 
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APPENDIX B: POISSON SUMMATION FORMULAS 



In this Appendix we recall the well-known Poisson summation formula and adapt it for functions with finite support 
in K. In the theory of generalized functions |3l| the following identity is proven: 

00 00 

J2 S{x-m)= e^™ (Bl) 

m— — 00 m— — CO 



Let / : [a,P] — > M be a continuous function with continuous derivative on the interval [q:,/3]. Multiply Eq. (Bl) 
on both sides with f{x), then integrate both sides from a to /?. In the evaluation of the left hand side we have to 
pay attention to whether any, or both the numbers a and (3 are integers, or non-integers. In the integer case the 
contribution of the end-point is calculated via the identity: 

n+r 

J dx6ix-n) fix) = ^f{n) , Vr>0 (B2) 

n 

Assuming that / is absolutely integrable if /3 = cx), and choosing a — 0, the classical Poisson summation formula is 
obtained: 

oo oo 

I] /(^) = 2'^^°) + / "^"^ + ^ ^ / "^"^ -^("^^ cos(2™) (B3) 
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Let us write also explicitely out the case when both a and (3 are integers: 



/('^) = 2^^^"^) + /(Z^)] + y /(^) + 2 E y /(^) cos(27rmx) , when a,/3 e Z 



(B4) 



Next we apply these equations to give an exact closed expression for the slope correlation function for finite L [eq. 
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We have \a\ < 1, and 
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In order to apply the Poisson summation formula (B4), we introduce the function: 
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and identify in (B4) a = 1 and /3 = L — 1. The non-integral terms of ( p34| ) give: 

1, 



2 2k-L 1 — a cos [^ij 



The next term becomes: 
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where during the evaluation of the integral we made a simple change of variables and used a well-known integral from 
random walk theory m , [p2| : 



7r 

Jdx- 



2tt 



— a cos X yT 



I > 



(Bll) 



The sum over the integrals in (B4) can also be evaluated, and one obtains: 



L-l 



2 J dx f{x) cos(27rna:) = 



a {U + 6-') 



2it/L 



2kVT^ 1-^^ 



> / da; cosfnia;) — 

2ttk^ J ^ ' 1 



Alx 



acosx 



-2n/L 



where 



and |6| < 1 



(B12) 



(B13) 



To compute the sum on the rhs of (B12) we recall another identity from the theory of generalized functions (see Ref. 
P, page 155): 



^e-=^ ^ 5{x~2m^^)+'-ctg[^) 



m— — OO 



x\ 1 
2 



(B14) 
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Combining (B14) and identity (Bl), one obtains: 

oo oo 

cos(na:) = tt d{x — 2m7r) + — 



m— — oo 



Peforming the sum over n directly in the rhs of (B12) via (pla), yields 



2tt/L 



a 

2HK 



E 



ilx 



dx 



co8{nLx)e 
1 — a cos X 



a X ^ 

77. 2^ 



2kL 



-2-!t/L 



m— — oo 



dy 



'^yS{y - 2TO7r) a 



Only m = ±1, contribute in ( |B16D . With the help of (|B2| ): 

2ir/L 



a ^-^ / , cos(nLx)e*'^ a 
^E / 



1 



^0 



2it/L 



1 — a cosy 27rK 



dx 



cos /x 
1 — a cos X 



2Tr/L 



1 — a cos a; 



-2ir/L 



2kL 1 1 - a 1 - acos 



27rK 



da; 



cos Za: 
1 — a cos a; 



(B15) 



(B16) 



(B17) 



Using ( [B17| ) in (|B12 ), we can add the result to the rest of the contributions (B£) and ( p310| ) to obtain the final 
expression [eq. (oM] after the cancellations. 
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